pro gauge_ch4_intercomparison_zonal_plot, z, x, y, filename, level, title

  vertPos=[0.15, 0.15, 0.8, 0.9]
  cbPos=[0.85, 0.1, 0.9, 0.9]

  cgLoadCT, 26, NColors=n_elements(level), Bottom=1, Brewer=1
  
  cgPS_Open, strcompress(filename + '.ps', /remove_all)
  
  ; Draw the contour plot, without Y axes, which we add later.
  cgContour, z, x, y, XTitle='Latitude', XStyle=1, $
    YTitle='Pressure (hPa)', YStyle=1, $
    Position=vertPos, Label=1, /Fill, /Outline, $
    C_Colors=Indgen(n_elements(level)), Levels=level, $
    yRange=[1000., 10.], title=title
  
  ; Add the colorbar.
  cgColorbar, Divisions=n_elements(level), NColors=n_elements(level)-1, Bottom=1, Range=[min(level), max(level)], $
    Title='Zonal mole fraction anomaly (ppb)', $
    TLocation='right', /vertical, /Discrete, Position=cbPos
    
  cgPS_Close, filetype='PNG', /delete_ps, density=200

end


pro gauge_ch4_intercomparison

  species='CH4'

  ;Get TOMCAT
  fname=file_search(strcompress('/data/shared/GAUGE/CH4/TOMCAT/CH4_N2O_TOMCAT3D_ppb_global_1985_2011.nc', /remove_all))

  lonTOMCAT=mr_ncget('lon', filename=fname)
  latTOMCAT=mr_ncget('lat', filename=fname)
  dateTOMCAT=mr_ncget('date', filename=fname)
  concTOMCAT=mr_ncget('mf_ch4', filename=fname)
  hyam=mr_ncget('hyam', filename=fname)
  hybm=mr_ncget('hybm', filename=fname)

  lonSize=n_elements(lonTOMCAT)
  latSize=n_elements(latTOMCAT)
  levSize=n_elements(hyam)
  timeSize=n_elements(dateTOMCAT)

  PTOMCAT=fltarr(LonSize, latSize, levSize, timeSize)

  PS=mr_ncget('PS', filename=fname)
  for ti=0, timeSize-1 do begin
    PTOMCAT[*, *, *, ti]=mr_hybrid_coords(hyam, hybm, PS[*, *, ti])
  endfor

  levTOMCAT=fltarr(levSize)
  for levI=0, levSize-1 do begin
    levTOMCAT[levI]=mean(PTOMCAT[*, *, levI, *])/100.
  endfor

  ;Get GEOS-Chem
  fnames=file_search(strcompress('/data/shared/GAUGE/CH4/GEOS-Chem/*.nc', /remove_all))
  fname=fnames[0]

  lonGeos=mr_ncget('longitude', filename=fname)
  latGeos=mr_ncget('latitude', filename=fname)
  hyam=mr_ncget('ab', filename=fname)
  hybm=mr_ncget('bp', filename=fname)

  lonSize=n_elements(lonGeos)
  latSize=n_elements(latGeos)
  levSize=n_elements(hyam)-1
  timeSize=n_elements(fnames)

  concgeos=fltarr(LonSize, latSize, levSize, timeSize)
  Pgeos=fltarr(LonSize, latSize, levSize, timeSize)
  dateGeos=lonarr(timeSize)
  
  foreach fname, fnames, ti do begin
    PS=mr_ncget('pres', filename=fname)
    Pgeos[*, *, *, ti]=reverse(transpose(mr_hybrid_coords(hyam, hybm, PS, /half), [1, 0, 2]), 3)
    concGeos[*, *, *, ti]=reverse(transpose(mr_ncget('Concentration', filename=fname), [2, 1, 0]), 3)*1.e9
    dateGeos[ti]=mr_ncget('date', filename=fname)
  endforeach

  levGeos=fltarr(levSize)
  for levI=0, levSize-1 do begin
    levGeos[levI]=mean(pGeos[*, *, levI, *])
  endfor

  
  ;Get MOZART
  mzOut=gauge_mozart_read_history('CH4', 'CH4', 'CH4_VMR_avrg', 1996, last_year=2010)
  mzOut.var=mzOut.var*1.e9

  locations=hash()
  locations['SPO']=[-89.98, -24.8]
  locations['MHD']=[53.326, -9.899]
  locations['MLO']=[19.5362, -155.5763]
  locations['SMO']=[-14.2474, -170.5644]

  keys=locations.keys()
  
  foreach key, keys do begin
    ;MOZART
    whLat=mr_nearest(mzOut.lat, locations[key, 0])
    whLon=mr_nearest(mzOut.lon, locations[key, 1])
    xMoz=mr_timedec(mr_mozdate(mzOut.date)-15.)
    yMoz=reform((mzOut.var)[whLon, whLat, -1, *])
    outMoz=mr_kz_filter(xMoz, yMoz, window=0.7, iterations=4, /growth)

    yS=mr_kz_filter(xMoz, yMoz, window=0.7, iterations=4)
    yDetrend=yMoz - interpol(yS.y, yS.x, xMoz)
    yseasonMoz=fltarr(12)
    counter=fltarr(12)
    
    ; Extract month for calculating the seasonal cycle
    moz_month=fix(12*(xMoz - fix(xMoz))) + 1
    
    for mi=0, 11 do begin
      month_index = where(Moz_month eq mi +1, count)
      if count gt 0 then yseasonMoz[mi]=mean(yDetrend[Month_index], /nan) else yseasonMoz[mi] = !values.F_NAN
    endfor

    ;GEOS-Chem
    whLat=mr_nearest(latGeos, locations[key, 0])
    whLon=mr_nearest(lonGeos, locations[key, 1])
    xGeos=mr_timedec(mr_mozdate(dateGeos))
    yGeos=reform((concGeos)[whLon, whLat, -1, *])
    outGeos=mr_kz_filter(xGeos, yGeos, window=0.7, iterations=4, /growth)

    yS=mr_kz_filter(xGeos, yGeos, window=0.7, iterations=4)
    yDetrend=yGeos - interpol(yS.y, yS.x, xGeos)
    yseasonGeos=fltarr(12)
    
    ; Extract month for calculating the seasonal cycle
    Geos_month=fix(12*(xGeos - fix(xGeos))) + 1

    for mi=0, 11 do begin
      month_index = where(Geos_month eq mi +1, count)
      if count gt 0 then yseasonGeos[mi]=mean(yDetrend[Month_index], /nan) else yseasonGeos[mi] = !values.F_NAN
    endfor

    ;TOMCAT
    whLat=mr_nearest(latTOMCAT, locations[key, 0])
    whLon=mr_nearest(lonTOMCAT, locations[key, 1])
    xTOMCAT=mr_timedec(mr_mozdate(dateTOMCAT))
    yTOMCAT=reform((concTOMCAT)[whLon, whLat, -1, *])
    outTOMCAT=mr_kz_filter(xTOMCAT, yTOMCAT, window=0.7, iterations=4, /growth)

    yS=mr_kz_filter(xTOMCAT, yTOMCAT, window=0.7, iterations=4)
    yDetrend=yTOMCAT - interpol(yS.y, yS.x, xTOMCAT)
    yseasonTOMCAT=fltarr(12)
    
    ; Extract month for calculating the seasonal cycle
    TOMCAT_month=fix(12*(xTOMCAT - fix(xTOMCAT))) + 1
    
    for mi=0, 11 do begin
      month_index = where(TOMCAT_month eq mi +1, count)
      if count gt 0 then yseasonTOMCAT[mi]=mean(yDetrend[Month_index], /nan) else yseasonTOMCAT[mi] = !values.F_NAN
    endfor

    ;MEASUREMENTS
    xNOAA=mr_ncget('time', $
      filename='/data/shared/GAUGE/CH4/NOAA_flask/processed/ch4_NOAA_CCGG_' + key + '_flask.nc', $
      attributes=xAtt)
    xREf=mr_time_convert(strmid(xAtt[where(xAtt[*, 0] eq 'units'), 1] , 14), 'YYYY-MM-DD HH:MM:SS')
    caldat, xRef, mNOAA, dNOAA, yNOAA
    xNOAA=julday(replicate(mNOAA, n_elements(xNOAA)), $
      replicate(dNOAA, n_elements(xNOAA)), $
      replicate(yNOAA, n_elements(xNOAA)), 0, 0, xNOAA)

    yNOAA=mr_ncget('mole_fraction_of_methane_in_air', filename='/data/shared/GAUGE/CH4/NOAA_flask/processed/ch4_NOAA_CCGG_' + key + '_flask.nc')
    combined=mr_combine_measurements(xNOAA, yNOAA, /monthly)
    xNOAA=mr_timedec(combined.T)
    yNOAA=combined.x
    wh=where(xNOAA gt float(fix(xNOAA[0])+1) and xNOAA lt float(fix(xNOAA[-1])))
    xNOAA=xNOAA[wh]
    yNOAA=yNOAA[wh]

    yS=mr_kz_filter(xNOAA, yNOAA, window=0.7, iterations=4)
    yDetrend=yNOAA - interpol(yS.y, yS.x, xNOAA)
    yseasonNOAA=fltarr(12)
    dyseasonNOAA=fltarr(12)
    
    ; Extract month for calculating the seasonal cycle
    NOAA_month=fix(12*(xNOAA - fix(xNOAA))) + 1

    
    for mi=0, 11 do begin
      month_index = where(NOAA_month eq mi +1, count)
      if count gt 0 then yseasonNOAA[mi]=mean(yDetrend[month_index], /nan) else yseasonNOAA[mi]=!values.F_NAN
      if count gt 0 then dyseasonNOAA[mi]=stddev(yDetrend[month_index], /nan) else dyseasonnNOAA[mi]=!values.F_NAN
    endfor

    ;Mole fraction plot
    pl=hash()
    pl[0]=plot(xNOAA, yNOAA, color='grey', symbol='o', thick=0.5, $
      ytitle='Mole fraction (ppb)', buffer=1, name='NOAA', linestyle=' ')
    pl[1]=plot(xMoz, yMoz, title=key, color='blue', thick=2, $
      position=[0.15, 0.5, 0.85, 0.95], name='MOZART', /over, /current)
    pl[2]=plot(xGeos, yGeos, color='red', thick=2, /over, /current, $
      ytitle='Mole fraction (ppb)', name='GEOS-Chem')
    pl[3]=plot(xTOMCAT, yTOMCAT, color='purple', thick=2, /over, /current, $
      ytitle='Mole fraction (ppb)', name='TOMCAT')
    leg=legend(/auto_text_color, target=pl, shadow=0, position=[0.4, 0.93])
    
    pl[0].xrange=[1988., 2012.]
    pl[0].yrange=[1600., 2200.]
    ax=pl[0].axes
    ax[0].hide=1
    
    ;Growth rate plot
    gpl=plot(outMoz.x, outMoz.y, color='blue', thick=2, $
      position=[0.15, 0.35, 0.85, 0.5], /current)
    gpl=plot(outGeos.x, outGeos.y, color='red', thick=2, /over, /current)
    gpl=plot(outTOMCAT.x, outTOMCAT.y, color='purple', thick=2, /over, /current)
    gpl[0].xrange=[1988., 2012.]
    ax=gpl.axes
    ax[1].hide=1
    ax[3].hide=1
    yAxis=axis('Y', location=['RIGHT'], target=gpl, $
      title='Growth rate (ppb/yr)')
    yAxis=axis('Y', location=['LEFT'], target=gpl, tickname=replicate(' ', 5))

    ;seasonal cycle rate plot
    spl=plot(indgen(12)+1, ySeasonMOZ, color='blue', thick=2, $
      position=[0.15, 0.1, 0.85, 0.25], /current, ytitle='Monthly !c mean (ppb)')
    spl=plot(indgen(12)+1, ySeasonGEOS, color='red', thick=2, /over, /current)
    spl=plot(indgen(12)+1, ySeasonTOMCAT, color='purple', thick=2, /over, /current)
    spl=errorplot(indgen(12)+1, ySeasonNOAA, dySeasonNOAA, color='grey', symbol='o', thick=0.5, $
      linestyle=' ', /over, /current)
    spl.xrange=[1, 12]
    spl.yminor=1

    ;Save plot
    pl[0].save, strcompress('/data/shared/GAUGE/CH4/Intercomparison/CH4_station_' + key + '.png', /remove_all), $
      resolution=200
    pl[0].close
    
  endforeach

  ;Plot inter-hemispheric gradient, map and vertical gradient
  IHGMoz=hash()
  IHGGeos=hash()
  IHGTOMCAT=hash()
  IHGMoz['jan']=fltarr(n_elements(mzOut.lat))
  IHGMoz['jul']=fltarr(n_elements(mzOut.lat))
  IHGGeos['jan']=fltarr(n_elements(latGeos))
  IHGGeos['jul']=fltarr(n_elements(latGeos))
  IHGTOMCAT['jan']=fltarr(n_elements(latTOMCAT))
  IHGTOMCAT['jul']=fltarr(n_elements(latTOMCAT))

  mapMoz=hash()
  mapGeos=hash()
  mapTOMCAT=hash()
  mapMoz['jan']=fltarr(n_elements(mzOut.lon), n_elements(mzOut.lat))
  mapMoz['jul']=fltarr(n_elements(mzOut.lon), n_elements(mzOut.lat))
  mapGeos['jan']=fltarr(n_elements(lonGeos), n_elements(latGeos))
  mapGeos['jul']=fltarr(n_elements(lonGeos), n_elements(latGeos))
  mapTOMCAT['jan']=fltarr(n_elements(lonTOMCAT), n_elements(latTOMCAT))
  mapTOMCAT['jul']=fltarr(n_elements(lonTOMCAT), n_elements(latTOMCAT))

  vertMoz=hash()
  vertGeos=hash()
  vertTOMCAT=hash()
  vertMoz['jan']=fltarr(n_elements(mzOut.lat), n_elements(mzOut.lev))
  vertMoz['jul']=fltarr(n_elements(mzOut.lat), n_elements(mzOut.lev))
  vertGeos['jan']=fltarr(n_elements(latGeos), n_elements(levGeos))
  vertGeos['jul']=fltarr(n_elements(latGeos), n_elements(levGeos))
  vertTOMCAT['jan']=fltarr(n_elements(latTOMCAT), n_elements(levTOMCAT))
  vertTOMCAT['jul']=fltarr(n_elements(latTOMCAT), n_elements(levTOMCAT))
  
  count=0
  AvStartY=2004
  AvEndY=2010
  nYears=AvEndY - AvStartY
  startYMoz=(where(mr_timedec(mr_mozdate(mzOut.date)-15.) ge float(AvStartY)))[0]
  startYGEOS=(where(mr_timedec(mr_mozdate(dateGeos)) ge float(AvStartY)))[0]
  startYTOMCAT=(where(mr_timedec(mr_mozdate(dateTOMCAT)) ge float(AvStartY)))[0]
  
  for yi=0, nYears-1 do begin
    whLon=mr_nearest(mzOut.lon, 180.)
    IHGMoz['jan']+=(mzOut.var)[whLon, *, -1, startYMoz + yi*12] - (mzOut.var)[whLon, 0, -1, startYMoz + yi*12]
    IHGMoz['jul']+=(mzOut.var)[whLon, *, -1, startYMoz + yi*12+6] - (mzOut.var)[whLon, 0, -1, startYMoz + yi*12+6]
    mapMoz['jan']+=(mzOut.var)[*, *, -1, startYMoz + yi*12] - (mzOut.var)[whLon, 0, -1, startYMoz + yi*12]
    mapMoz['jul']+=(mzOut.var)[*, *, -1, startYMoz + yi*12+6] - (mzOut.var)[whLon, 0, -1, startYMoz + yi*12+6]
    vertMoz['jan']+=mean((mzOut.var)[*, *, *, startYMoz + yi*12], dimension=1) - (mzOut.var)[whLon, 0, -1, startYMoz + yi*12]
    vertMoz['jul']+=mean((mzOut.var)[*, *, *, startYMoz + yi*12+6], dimension=1) - (mzOut.var)[whLon, 0, -1, startYMoz + yi*12+6]

    whLon=mr_nearest(lonGeos, 180.)
    IHGGeos['jan']+=(concGeos)[whLon, *, -1, startYGEOS + yi*12] - (concGeos)[whLon, 0, -1, startYGEOS + yi*12]
    IHGGeos['jul']+=(concGeos)[whLon, *, -1, startYGEOS + yi*12+6] - (concGeos)[whLon, 0, -1, startYGEOS + yi*12+6]
    mapGeos['jan']+=(concGeos)[*, *, -1, startYGEOS + yi*12] - (concGeos)[whLon, 0, -1, startYGEOS + yi*12]
    mapGeos['jul']+=(concGeos)[*, *, -1, startYGEOS + yi*12+6] - (concGeos)[whLon, 0, -1, startYGEOS + yi*12+6]
    vertGeos['jan']+=mean((concGeos)[*, *, *, startYGEOS + yi*12], dimension=1) - (concGeos)[whLon, 0, -1, startYGEOS + yi*12]
    vertGeos['jul']+=mean((concGeos)[*, *, *, startYGEOS + yi*12+6], dimension=1) - (concGeos)[whLon, 0, -1, startYGEOS + yi*12+6]

    whLon=mr_nearest(lonTOMCAT, 180.)
    IHGTOMCAT['jan']+=(concTOMCAT)[whLon, *, -1, startYTOMCAT+ yi*12] - (concTOMCAT)[whLon, 0, -1, startYTOMCAT+ yi*12]
    IHGTOMCAT['jul']+=(concTOMCAT)[whLon, *, -1, startYTOMCAT+ yi*12+6] - (concTOMCAT)[whLon, 0, -1, startYTOMCAT+ yi*12+6]
    mapTOMCAT['jan']+=(concTOMCAT)[*, *, -1, startYTOMCAT+ yi*12] - (concTOMCAT)[whLon, 0, -1, startYTOMCAT+ yi*12]
    mapTOMCAT['jul']+=(concTOMCAT)[*, *, -1, startYTOMCAT+ yi*12+6] - (concTOMCAT)[whLon, 0, -1, startYTOMCAT+ yi*12+6]
    vertTOMCAT['jan']+=mean((concTOMCAT)[*, *, *, startYTOMCAT+ yi*12], dimension=1) - (concTOMCAT)[whLon, 0, -1, startYTOMCAT+ yi*12]
    vertTOMCAT['jul']+=mean((concTOMCAT)[*, *, *, startYTOMCAT+ yi*12+6], dimension=1) - (concTOMCAT)[whLon, 0, -1, startYTOMCAT+ yi*12+6]

    count++
  endfor

  foreach key, ['jan', 'jul'] do begin
    IHGMoz[key]=IHGMoz[key]/float(count)
    IHGGeos[key]=IHGGeos[key]/float(count)
    IHGTOMCAT[key]=IHGTOMCAT[key]/float(count)
    mapMoz[key]=mapMoz[key]/float(count)
    mapGeos[key]=mapGeos[key]/float(count)
    mapTOMCAT[key]=mapTOMCAT[key]/float(count)
    vertMoz[key]=vertMoz[key]/float(count)
    vertGeos[key]=vertGeos[key]/float(count)
    vertTOMCAT[key]=vertTOMCAT[key]/float(count)
  endforeach

  ;Plot surface inter-hemispheric gradient
  pl=hash()
  pl[0]=plot(mzOut.lat, IHGMoz['jan'], color='blue', name='MOZART (January)', $
    ytitle='Latitudinal gradient (ppb)', /buffer)
  pl[1]=plot(mzOut.lat, IHGMoz['jul'], color='blue', linestyle='--', /current, /over, $
    name='MOZART (July)')
  pl[2]=plot(latGeos, IHGGeos['jan'], color='red', name='GEOS-Chem (January)', $
    /current, /over)
  pl[3]=plot(latGeos, IHGGeos['jul'], color='red', linestyle='--', /current, /over, $
    name='GEOS-Chem (July)')
  pl[4]=plot(latTOMCAT, IHGTOMCAT['jan'], color='purple', name='TOMCAT (January)', $
    /current, /over)
  pl[5]=plot(latTOMCAT, IHGTOMCAT['jul'], color='purple', linestyle='--', /current, /over, $
    name='TOMCAT (July)')
  leg=legend(/auto_text_color, target=pl, shadow=0, position=[0.5, 0.8])
  lpl=plot([-100., 100.], [0., 0.], /current, /over)
  pl[0].save, strcompress('/data/shared/GAUGE/CH4/Intercomparison/CH4_lat_gradient.png')
  pl[0].close

  ;Plot map
  mapPos=[0.1, 0.1, 0.8, 0.9]
  cbPos=[0.85, 0.1, 0.9, 0.9]
  level=findgen(21)*20. - 100*2.
  RGB_indices=indgen(21)*255/20

  plmap=mapMoz['jan']
  x=mzOut.lon
  y=mzOut.lat
  mr_grid_orient, x, y, plmap, /neglon
  mpl=contour(plmap, x, y, /fill, RGB_table=73, $
    xstyle=1, ystyle=1, title='MOZART January', position=mapPos, /buffer, $
    c_value=level, rgb_indices=rgb_indices)
  c = COLORBAR(TARGET=mpl, ORIENTATION=1, $
    POSITION=cbPos, $
    TITLE='Mole fraction (ppb)', textpos=1)
  mpl.save, '/data/shared/GAUGE/CH4/Intercomparison/CH4_map_MOZART_Jan.png', resolution=200
  mpl.close

  plmap=mapMoz['jul']
  x=mzOut.lon
  y=mzOut.lat
  mr_grid_orient, x, y, plmap, /neglon
  mpl=contour(plmap, x, y, /fill, RGB_table=73, $
    xstyle=1, ystyle=1, title='MOZART July', position=mapPos, /buffer, $
    c_value=level, rgb_indices=rgb_indices)
  c = COLORBAR(TARGET=mpl, ORIENTATION=1, $
    POSITION=cbPos, $
    TITLE='Mole fraction (ppb)', textpos=1)
  mpl.save, '/data/shared/GAUGE/CH4/Intercomparison/CH4_map_MOZART_Jul.png', resolution=200
  mpl.close

  plmap=mapGeos['jan']
  x=lonGeos
  y=latGeos
  y[0]=-90
  y[-1]=90
  mpl=contour(plmap, x, y, /fill, RGB_table=73, $
    xstyle=1, ystyle=1, title='GEOS-Chem January', position=mapPos, /buffer, $
    c_value=level, rgb_indices=rgb_indices)
  c = COLORBAR(TARGET=mpl, ORIENTATION=1, $
    POSITION=cbPos, $
    TITLE='Mole fraction (ppb)', textpos=1)
  mpl.save, '/data/shared/GAUGE/CH4/Intercomparison/CH4_map_GEOS-Chem_Jan.png', resolution=200
  mpl.close

  plmap=mapGeos['jul']
  x=lonGeos
  y=latGeos
  y[0]=-90
  y[-1]=90
  mpl=contour(plmap, x, y, /fill, RGB_table=73, $
    xstyle=1, ystyle=1, title='GEOS-Chem July', position=mapPos, /buffer, $
    c_value=level, rgb_indices=rgb_indices)
  c = COLORBAR(TARGET=mpl, ORIENTATION=1, $
    POSITION=cbPos, $
    TITLE='Mole fraction (ppb)', textpos=1)
  mpl.save, '/data/shared/GAUGE/CH4/Intercomparison/CH4_map_GEOS-Chem_Jul.png', resolution=200
  mpl.close

  plmap=mapTOMCAT['jan']
  x=lonTOMCAT
  y=latTOMCAT
  y[0]=-90
  y[-1]=90
  mpl=contour(plmap, x, y, /fill, RGB_table=73, $
    xstyle=1, ystyle=1, title='TOMCAT January', position=mapPos, /buffer, $
    c_value=level, rgb_indices=rgb_indices)
  c = COLORBAR(TARGET=mpl, ORIENTATION=1, $
    POSITION=cbPos, $
    TITLE='Mole fraction (ppb)', textpos=1)
  mpl.save, '/data/shared/GAUGE/CH4/Intercomparison/CH4_map_TOMCAT_Jan.png', resolution=200
  mpl.close

  plmap=mapTOMCAT['jul']
  x=lonTOMCAT
  y=latTOMCAT
  y[0]=-90
  y[-1]=90
  mpl=contour(plmap, x, y, /fill, RGB_table=73, $
    xstyle=1, ystyle=1, title='TOMCAT July', position=mapPos, /buffer, $
    c_value=level, rgb_indices=rgb_indices)
  c = COLORBAR(TARGET=mpl, ORIENTATION=1, $
    POSITION=cbPos, $
    TITLE='Mole fraction (ppb)', textpos=1)
  mpl.save, '/data/shared/GAUGE/CH4/Intercomparison/CH4_map_TOMCAT_Jul.png', resolution=200
  mpl.close



  ;Plot zonal profiles
  level=findgen(21)*20. - 200.

  gauge_CH4_intercomparison_zonal_plot, vertMoz['jan'], mzOut.lat, mzOut.lev, $
    '/data/shared/GAUGE/CH4/Intercomparison/CH4_vert_MOZART_jan', level, $
    'MOZART, (Jan, 180E, difference from SPO)'
  gauge_CH4_intercomparison_zonal_plot, vertMoz['jul'], mzOut.lat, mzOut.lev, $
    '/data/shared/GAUGE/CH4/Intercomparison/CH4_vert_MOZART_jul', level, $
    'MOZART, (July, 180E, difference from SPO)'
  gauge_CH4_intercomparison_zonal_plot, vertGeos['jan'], latGeos, levGeos, $
    '/data/shared/GAUGE/CH4/Intercomparison/CH4_vert_GEOS-Chem_jan', level, $
    'GEOS-Chem, (Jan, 180E, difference from SPO)'
  gauge_CH4_intercomparison_zonal_plot, vertGeos['jul'], latGeos, levGeos, $
    '/data/shared/GAUGE/CH4/Intercomparison/CH4_vert_GEOS-Chem_jul', level, $
    'GEOS-Chem, (July, 180E, difference from SPO)'
  gauge_CH4_intercomparison_zonal_plot, vertTOMCAT['jan'], latTOMCAT, levTOMCAT, $
    '/data/shared/GAUGE/CH4/Intercomparison/CH4_vert_TOMCAT_jan', level, $
    'TOMCAT, (Jan, 180E, difference from SPO)'
  gauge_CH4_intercomparison_zonal_plot, vertTOMCAT['jul'], latTOMCAT, levTOMCAT, $
    '/data/shared/GAUGE/CH4/Intercomparison/CH4_vert_TOMCAT_jul', level, $
    'TOMCAT, (July, 180E, difference from SPO)'
  
  spawn, '/data/shared/GAUGE/CH4/Intercomparison/combine_plots'
  
end